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1.  Summary 


In  this  document  we  report  on  activities  related  to  the  project  “Anisoplanatic  Imaging  through 
Turbulence”,  grant  number  FA9550-14-1-0244.  The  activities  in  Year  2  have  followed  two  main  research 
directions: 

The  development  of  a  new  point-spread  function  (PSF)  reconstruction  approach  which  aims  at 
extracting,  from  multi-frame  observations,  a  number  of  PSFs  at  several  locations  within  an  image 
of  an  arbitrary  object  (without  point  sources)  (Section  3.1). 

The  development  of  a  marginal  blind  deconvolution  estimator,  combined  with  constraints  on  the 
observed  object  such  as  positivity,  for  identification  of  long-exposure  adaptive  optics  PSFs.  This 
approach  can  also  be  extended  to  the  anisoplanatic  case  by  its  application  in  different  locations  of 
the  field  of  view  as  a  first  building  block  of  a  global  shift-variant  image  restoration  approach 
(Section  3.2). 

Both  approaches  are  complementary  to  each  other,  e.g.,  from  the  first  approach  it  is  possible  to  extract  a 
model  for  the  object’s  power  spectral  density  that  can  be  used  to  complete  the  parameters  to  be  used  in  the 
second  approach. 

2.  Introduction 

The  performance  of  optical  systems  is  degraded  by  atmospheric  turbulence  when  observing  vertically  (e.g. 
astronomy)  or  horizontally  (e.g.  surveillance,  military  reconnaissance).  This  degradation  can  be  alleviated 
in  software  (real-time  de -blurring  or  post-processing)  or  hardware  (adaptive  optics  -  AO).  Tremendous 
progress  has  been  achieved  in  this  area:  in  astronomy  almost  every  major  observatory  is  now  equipped 
with  first-generation  AO  systems  and  some  second-generation  systems  are  currently  coming  online.  High- 
resolution  imaging  of  satellites  is  now  a  routine  task  at  the  Air  Force  Maui  Optical  and  Supercomputing 
(AMOS)  observatory  thanks  to  a  combination  of  AO  and  post -processing.  The  technology  and  associated 
image  processing  techniques  have  been  successfully  transferred  to  commercial  markets,  especially  in  the 
field  of  ophthalmic  imaging. 

One  problem  still  remains:  Adaptive  optics  systems  are  generally  only  capable  of  correcting  very  small 
fields  of  view.  In  single-conjugated  AO  systems,  i.e.,  ones  that  measure  and  correct  atmospheric 
distortions  in  one  direction,  with  a  single  deformable  mirror,  the  usable  field-of-view  is  determined  by  the 
so-called  isoplanatic  angle  (Figure  1).  This  angle  describes  the  field-of-view  in  which  turbulence-induced 
aberrations  could  be  considered  constant  and  therefore  correctable  with  a  single  device.  In  other  words, 
given  a  certain  correction  direction,  isoplanatic  angle  gives  the  maximum  angular  separation  from  this 
direction  at  which  reasonably  good  correction  can  be  expected.  In  conventional  AO  systems,  the 
correctable  field-of-view  has  a  size  of  the  order  of  isoplanatic  angle. 

In  theory,  the  problem  of  wide-field  imaging  can  be  solved  by  multi-conjugate  AO.  Such  systems,  already 
deployed  successfully  at  a  handful  of  astronomical  observatories,  correct  aberrations  resulting  from 
propagations  over  distinct  paths  using  a  number  of  wavefront  sensors  and  corresponding  deformable 
mirrors.  These  systems  are  complex,  bulky,  one-off  instruments.  In  astronomical  imaging,  two-three 
wavefront  sensors  and  deformable  mirrors  already  increase  the  field  of  view  to  enable  new  science. 
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Reference  Star 


Target 


Z  altitude 


Z  altitude 


Fig.  1.  Illustration  of  the  problem  of  anisoplanatism  in  imaging  through  turbulence  with  adaptive  optics, 
vertical-path  scenario.  Left:  when  strong  atmospheric  layer  is  located  close  to  the  telescope  pupil  then 
there  exists  a  large  overlap  in  significant  wavefront  aberrations  between  light  coming  from  the  object  of 
interest  and  the  guide  star  (reference  for  AO).  Right:  when  the  layer  of  significant  turbulence  is  located 
further  up  then  the  overlap  between  the  wavefronts  coming  from  the  two  directions  is  reduced  and 
subsequently  the  correction  applied  to  the  deformable  mirror  is  valid  only  for  the  reference  star  and  not 
for  the  target. 

This  solution  is  not  applicable  for  at  least  two  areas  of  interest:  For  ground-level  surveillance,  where  the 
field  of  view  could  extend  over  hundreds  of  isoplanatic  patches  because  of  the  strength  of  turbulence 
along  horizontal  paths  near  the  ground,  multi-conjugate  AO  would  have  to  consist  of  hundreds  of 
wavefront  sensors  and  corresponding  deformable  mirrors.  This  poses  an  enormous  technical  problem  and 
the  idea  of  multi-conjugate  AO  in  this  context  is  currently  not  being  followed.  Also,  in  air-to-ground, 
“looking-down”  scenario,  the  ratio  of  sensor  aperture  to  the  coherence  length  of  the  atmosphere  is  around 
unity  (as  opposed  to  tens  or  hundreds  in  imaging  from  the  ground)  but  isoplanatic  angles  are  usually  very 
small  because  optical  aberrations  are  located  far  away  from  the  imaging  pupil. 

The  hardware  solution  to  the  problem  of  anisoplanatism,  i.e.  multi -conjugate  AO,  is  very  expensive. 
Therefore,  in  this  project  we  follow  two  software  approaches. 


Fig.  2.  Simulated  data  showing  a  stellar  field  imaged  with  single-conjugate  AO  (for  more  details  on  the 
simulation  see  [1]).  Left:  original  data,  false  color  and  stretched  scale.  Right:  result  of  image 
reconstruction  with  a  single  PSF  corresponding  to  the  guide  star,  ignoring  anisoplanatism.  It  can  be  seen 
how  image  quality  degrades  when  moving  away  from  the  guide  star. 


Data  from  single-conjugate  AO  systems  are  difficult  to  post-process  (deconvolve)  because  of  the  spatially  - 
variant  blur  (Figure  2).  This  means  that  each  point  in  an  image  is  blurred  by  a  slightly  different  point  - 
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spread  function  (PSF).  On  the  other  hand,  almost  all  deconvolution  algorithms  rely  on  PSF  being  constant 
in  the  field  of  view.  This  is  the  shift -invariant  paradigm  which  does  not  lead  to  optimal  results  (cf.  right 
panel  of  Figure  2).  In  the  following  sections  we  describe  two  solutions  which  we  followed  to  tackle  this 
problem. 

3.  Methods,  Assumptions,  and  Procedures 
3.1  PSF  reconstruction  directly  from  target  data 

In  the  first  year  of  the  project  work  was  performed  on  reformulation  of  imaging  problem  from  traditional 
shift-invariant  into  the  shift -variant  case,  development  of  simulations  of  anisoplanatic  imagery  and 
anisoplanatic  image  reconstruction  algorithms.  Gains  have  been  made  on  all  three  fronts  but  the 
underlying  assumption  was  that  the  algorithms  have  access  to  knowledge  about  the  form  of  the  PSF  at 
every  point  in  an  image.  In  the  second  year,  efforts  were  devoted  to  finding  an  approach  that  could  reveal 
these  PSFs  even  for  objects  without  features  such  as  edges  or  bright  points.  As  mentioned  in  the  grant 
proposal,  building  on  our  previous  work  on  object-cancelling  transformations  [2]  we  hoped  to  find  new 
transformations  which  would  allow  to  reduce  the  problem  to  two  equations  with  two  unknowns,  giving 
directly  the  PSF.  Ultimately,  this  proved  possible  and  good  results  have  been  obtained  (see  Section  4.1). 

The  imaging  equation  for  a  true  object  o(v)  convolved  with  a  PSF  h(v) ,  giving  the  recorded  image  i(v) 
is: 


i(x)  =  o(3c)  ®  h(v)  ( 1 ) 

where  0  denotes  the  convolution  operator  and  the  symbol  x  corresponds  to  two-dimensional  focal-plane 
coordinate.  Equation  (1),  transformed  into  Fourier  domain,  squared  and  averaged  over  M  images, 
becomes: 


(lI(")f}„  =|0(")f(|H(«)|2)i) 


(2) 


If  we  gain  access  to  the  AO  speckle  transfer  function  ^|H(w)|y  then  the  AO  modulation  transfer  function 
(MTF)  |(H  (m))|  can  be  directly  obtained  (see  Equation  (7)).  A  PSF  is  then  derived  by  Fourier  transforming 

the  MTF.  Although  the  phase  of  the  optical  transfer  function  (OTF)  (h(m))  is  lost  in  the  process,  we  show 
that  the  approach  still  yields  a  useful  PSF. 

In  order  to  break  the  symmetry  of  Equation  (2)  we  cancel  out  the  object  by  transforming  it  into  equations 
involving  moments  of|H(w)|2.  One  such  transformation,  termed  “Fourier  contrast”,  was  already 

discovered  by  us  before  [2].  It  is  given  by  standard  deviation  of  the  power  spectrum  computed  for  each 
frequency  separately,  divided  by  the  mean  value  of  the  power  spectrum  at  this  frequency: 


c(m2) 


|0(«)|2  ]j(l 

H(m) 

f)' 

|0(«)|2  ( 

H(n) 

•> 

(3) 
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If  the  object  is  static  during  the  observations  then  its  power  spectrum  cancels  out  and  we  are  left  only  with 
the  part  involving  random  aberrations.  Note  that  all  static  terms  disappear  from  the  result,  i.e.  also 
telescope  OTF  will  be  removed.  Therefore  in  order  to  obtain  realistic  PSF  one  has  to  multiply  the  result  of 
the  method  again  by  the  telescope  OTF,  preferably  one  which  includes  also  static  errors  of  the  optics. 
Some  practical  issues  related  to  this  process  are  discussed  in  Section  4.1. 

Naturally,  from  the  statistical  moments  of  the  AO  speckle  transfer  function  it  is  not  possible  to  derive  its 
mean  value.  The  second  transformation  that  cancels  out  the  static  object  is  Fourier  skewness  [3]: 


■s(K«r) 


jocgjT  (lHf 

|0(m)|6 


+  2  H 


(4) 


The  transformation  amounts  to  computing  the  third  standardized  moment  of  the  image  power  spectra, 
separately  for  each  frequency.  If  the  object  does  not  change  during  the  observations  then  it  will  cancel  out 
from  Equation  (4). 

Equations  (3)  and  (4)  involve  many  higher-order  moments  of\H(u)\\  In  order  to  arrive  at  a  tractable 

formulation  of  the  problem  we  must  express  these  moments  using  as  few  parameters  as  possible.  Although 
in  general,  the  AO  OTF  is  modeled  with  many  variables  corresponding  to  the  parameters  of  the  turbulence 
and  the  AO  system  [4-6],  here  we  use  a  much  simpler  model.  We  start  with  the  realization  that  the  OTF  is 
a  speckle  [2]: 


,r>i  +A  ZUx  ,rii  +A  ZUy^  (5) 

where  <^,77  are  the  pupil-plane  coordinates,  A  is  the  average  wavelength  of  the  observations,  zis  the 
distance  from  the  exit  pupil  to  the  image  plane,  N  is  the  number  of  OTF  cells  in  the  area  of  overlap  [2]  and 
6  ’s  are  the  phases  -  identically  distributed  random  variables. 

Note  that  because  phases  are  random  the  OTF  also  has  to  be  random,  even  for  long  exposures.  Given  this 
realization,  one  can  write  moments  of  |H(w)|  using  only  the  following  parameters:  A,  which  is  of  little 

practical  interest,  and  various  powers  of  the  OTF  [7].  These  expressions  are  very  long  but  fortunately 
approximations  have  been  derived  for  the  case  of  large  N  [8].  For  contrast  and  skewness  one  obtains, 
respectively: 


H(m)  = 


s  « 


N»l 


(6) 
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Looking  at  Equation  (6)  we  see  that  we  obtained  two  equations  with  two  unknowns:  the  modulation 
transfer  function  |H|  and  the  number  of  OTF  cells  N.  One  can  now  proceed  to  calculate  modulation 
transfer  function  for  any  object  being  imaged.  Beforehand,  one  last  simplification  in  Equation  (6)  is  used: 


(fB-iuf)  =  [l  +  (AT  - l)exp(- De(A  zux,A  z%))]~exp(- D0(A  zux,A  zuY))&  |(h(m))|  ,  N  »1  (7) 

where D0() is  the  structure  function  of  the  AO  residual  phase.  This  approximation  is  only  valid  for  AO- 
compensated  imaging. 

To  summarize,  the  new  approach  to  AO  PSF  reconstruction  involves  the  recording  of  several  exposures  of 
the  same  object,  computing  contrast  and  skewness  of  the  power  spectra  of  the  data,  obtaining  the  AO  MTF 
from  Equations  (6)  and  (7)  and  Fourier-transforming  it  to  derive  the  average  PSF  of  the  observations.  As 
mentioned  before,  unfortunately  the  OTF  phase  information  is  lost  in  the  process. 


3.2  Marginal  blind  deconvolution  for  satellite  identification 

Blind  deconvolution,  both  in  the  astronomical  context  and  for  satellite  identification,  is  usually  performed 
in  a  joint  fashion,  i.e.,  object  and  PSF  are  estimated  together  alternating  the  minimization  of  a  certain  cost 
function  between  one  and  the  other.  However,  it  is  now  well  known  that  this  often  produces  a  solution  that 
is  degenerate  and,  when  it  works,  it  is  thanks  to  some  constraints  and  prior  information  about  the  PSF 
and/or  the  object,  such  as  positivity,  which  help  algorithms  converge  to  a  reasonable  solution. 

In  a  nutshell,  joint  blind  deconvolution  consists  in  maximizing  the  joint  probability  p(i,o,h),  where  i  is  the 
image  or  data,  and  o  and  h  are  the  unknown  object  and  PSF  to  be  recovered.  If  the  PSF  is  considered  a 
linear  combination  of  a  set  of  modes  {hj  }j=i ...j,  i.e.,  h=  T/j=iaj^]  >  where  the  set  of  {a}  are  the 
coefficients  that  parameterize  the  PSF,  then  the  joint  probability  to  be  maximized  would  be  p(i,o,a).  This 
change  of  parameter  (from  h  to  a)  allows  us  to  reduce  the  number  of  unknowns  from  the  number  of  pixels 
on  the  PSF  support  to  a  few  coefficients.  If  the  modes  are  also  PSFs  then  the  set  of  coefficients  is 
normalized  to  1  (Z]=±  aj  —  1)  and  each  of  them  is  positive  ( (Xj  >  0). 

In  the  case  where  the  noise  is  assumed  stationary  white  Gaussian,  and  a  Gaussian  prior  distribution  with 
mean  value  om  is  also  assumed  for  the  object,  then  the  joint  criterion  to  be  minimized  Jj(0mi)MAP>  which  is 
the  opposite  of  the  logarithm  of  the  joint  probability  p(i,o,a),  adopts  the  following  expression  [9]: 


J jMAP  C°f  a) 


ly  \I(u)-H(u)Q(u)\2  ly  \Q(u)-Om(u)\2 
2LU  Sn  ^  2LU  S0& ) 


(8) 


where  u  is  a  two-dimensional  spatial  frequency  coordinate,  Sn  and  S0  are,  respectively,  the  noise  and 
object  power  spectral  densities  (PSD),  and  capital  letters  denote  the  2-D  Fourier  Transform.  In  Equation 
(8)  if  we  replace  the  object  by  the  Wiener  solution  for  a  given  set  of  {a}  we  obtain  a  new  expression  that 
does  not  depend  explicitly  on  the  object. 


JjMAp(°w(a)>a ) 


ly  1  \I(u)-H(u)Om(u)\2 

~2^s0m  i hgop+s^ 


(9) 
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It  is  important  to  notice  that  JjMAp(°> a )  and  JjMAp(°w(a)> are  equivalent,  i.e.,  the  use  of  the  Wiener 
solution  does  not  change  the  shape  of  the  criterion  but  only  the  speed  at  which  it  is  minimized.  From 
Equation  (9)  it  is  easy  to  understand  why  the  joint  solution  is  degenerate:  if  the  mean  object  om  is 
constant,  and  since  the  PSF  as  well  as  the  set  of  coefficients  {a}  are  normalized  to  1,  the  numerator  does 
not  depend  on  such  set  of  coefficients  {a}.  Minimizing  J jMAp  is  equivalent  to  maximizing  the 
denominator,  i.e.,  choosing  the  PSF  with  the  highest \H(u)\,  which  is  the  sharpest  PSF.  Therefore  the 
solution  of  the  optimization  of  this  joint  criterion  is  always  the  sharpest  PSF,  whatever  the  true  PSF  may 
be! 

A  possible  solution  to  this  problem  is  the  use  of  the  marginal  estimator.  Its  goal  is  to  estimate  the  set  of 
most  likely  coefficients  {a}  on  average  for  all  possible  objects  by  integrating  the  object  o  out  of  the 
problem,  i.e.,  we  integrate  the  joint  probability  on  all  possible  objects.  In  other  words,  the  joint  estimation 
consists  in  maximizing  p(i,o,a),  whereas  the  marginal  estimation  consists  in  maximizing  j0p(i,o,a),  i.e., 
maximizing  the  joint  probability  over  all  possible  objects,  which  is  the  marginal  probability  p(i,a),. 
Marginalization  reduces  the  number  of  unknowns  providing  an  expression  that  only  depends  on  the  set  of 
PSF  coefficients  {a}.  Once  these  coefficients  are  estimated  they  are  used  to  restore  the  object,  for  instance 
by  Wiener  filtering  of  the  image.  If  we  again  assume  a  Gaussian  distribution  for  noise,  and  no  prior 
information  for  the  set  of  {a},  then  the  expression  of  the  marginal  criterion  is  [9]: 


7ml  (a)  =  2  Si 


1  \nu)-H(u)Om(u)\ 

S0  (v)  \H(u)\2  + 


^  +  i£aln(|H(S)l2+^)  +  C 

S0(u) 


(10) 


Equation  (10)  is  very  similar  to  Equation  (9)  however  the  presence  of  \H(u)\2  in  the  additional  term 
provides  the  marginal  criterion  with  better  properties  with  respect  to  the  joint  one,  as  can  be  seen  in  Figure  3. 


Alpha=0.30  Alpha-' 


Joint  criterion  for  0  <  a  <  1  Marginal  criterion  for  0  <  a  <  1 

Fig.  3.  Plot  for  the  joint  (left)  and  the  marginal  (right)  estimator  when  the  PSF  is  a  linear  combination  of 
two  modes  related  with  one  single  coefficient  of  value  a=0.3,  i.e.,  h  =  a  ■  hf0C  +  (1  —  a)  ■  hdef0C.  The 
reader  can  notice  that  the  marginal  criterion  has  its  minimum  at  the  correct  value  of  a,  whereas  the  joint 
criterion  has  its  minimum  for  that  a  that  pushes  the  PSF  to  be  the  focused  one,  i.e.,  the  sharpest  one. 
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From  Equations  (8)  and  (10)  is  easy  to  realize  that  the  marginal  and  the  joint  criteria  are  related  when  the 
Wiener  solution  is  used  for  the  latter: 


/ml  (a)  =  ]jMAp(.ow(d),a)  +  (|//(u)|2  +  C  (11) 

If  we  do  not  impose  the  Wiener  solution  in  the  joint  criterion,  we  can  create  a  new  criterion  derived  from 
the  marginal  one,  but  with  the  explicit  presence  of  the  object: 

Jml(°>  a)  =JjMAp(o,a)  +  (|//(w)|2  +  C  (12) 

This  new  criterion  can  be  minimized  alternating  between  the  object  and  the  PSF  coefficients,  and  the 
explicit  presence  of  the  object  allows  us  to  impose  strong  constraints  over  it,  such  as  positivity,  during  the 
optimization.  In  practice,  the  positivity  constraint  behaves  as  an  object  support  constraint  and  it  is 
especially  useful  in  the  astronomical  context  and  for  satellite  identification,  where  objects  to  be  recovered 
usually  lie  over  a  black  background.  The  drawback  of  using  the  non-analytical  solution  in  JjMAp(°>a) 
instead  of  the  Wiener  one  is  a  decrease  in  the  speed  of  convergence,  since  the  minimization  is  again 
performed  alternately  between  the  object  and  the  PSF. 

Then,  the  final  expression  for  the  new  marginal  criterion  with  positivity  on  the  object  is: 


Jml(o,  a)|s.t  o>0  2^^ 


\I(u)-H(u)0(u)\2 


|Q(g)-Om(u)|2 

S0(u) 


+ 


+  C  (13) 


Where  the  first  term  in  Equation  (13)  is  the  data  fidelity  term  when  stationary  white  Gaussian  noise  is 
considered,  the  second  term  is  a  Gaussian  prior  distribution  for  the  object,  with  an  explicit  presence  of  it, 
and  the  third  term  is  due  to  the  marginalization. 

Finally,  the  marginal  criterion,  either  with  or  without  positivity,  needs  the  object  and  noise  PSDs  S0(u ) 
and  Sn.  Fortunately,  both  can  be  estimated  together  with  the  PSF  coefficients  in  an  automatic  manner 
during  the  minimization.  In  order  to  reduce  the  number  of  new  parameters  to  estimate,  a  simple  model  for 

/c 

the  object  PSD  was  chosen,  e.g.,  S0  = - —p.  This  model  has  been  previously  used  with  success  in  both 

i+m 

\u0J 

the  astronomical  context  and  in  ground-based  imaging  of  satellites  [10],  and  it  only  adds  3  new  parameters 
{k,n0,p}  together  with  the  noise  PSD,  which  is  constant  since  we  are  assuming  white  noise. 

4.  Results  and  Discussion 

4.1  PSF  reconstruction  applied  to  simulated  and  real  data 

In  order  to  test  the  PSF  reconstruction  method  described  in  Section  3.1  we  first  performed  tests  on  images 
where  the  ground-truth  average  PSF  was  known  with  a  high  degree  of  accuracy:  a  large  set  of  AO  PSFs 
with  high  signal-to-noise  ratio.  A  sequence  of  4000  short  exposure  images  (texp  =  22  ms)  of  a  single  bright 
star  taken  with  the  3-m  Shane  Telescope  at  the  Lick  Observatory  was  used.  Images  were  taken  in  K-band 
(2.2  pm)  with  AO  (7x7  actuators  on  the  deformable  mirror)  switched  on.  The  resulting  Strehl  ratio  was 
53%.  Power  spectra  were  computed  for  each  image  and  subsequently  contrast  and  skewness  values  of 
each  frequency  component  were  computed,  i.e.  statistical  quantities  were  averaged  over  time. 
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As  mentioned  in  Section  3.1,  the  proposed  operations  on  the  data  remove  all  static  components  of  the 
images,  including  the  ones  that  must  be  modelled  accurately  to  capture  image  degradation.  While  it  is 
straightforward  to  re-introduce  the  OTF  of  the  telescope  with  an  arbitrarily  shaped  pupil  back  into  the 
result  of  the  method,  it  is  not  that  straightforward  to  recover  the  information  about  the  static  aberrations 
that  gets  lost  in  the  process.  This  is  an  issue  common  to  all  PSF  reconstruction  approaches  [6].  Here  we 
use  a  PSF  obtained  from  fiber  measurements  with  the  telescope  dome  closed,  i.e.  a  PSF  only  affected  by 
static  optical  aberrations.  This  PSF  is  convolved  with  the  AO-residual  PSF,  or  alternatively  their 
respective  OTFs  are  multiplied. 

The  azimuthally-averaged  results  of  applying  the  described  method  are  shown  in  the  left  panel  of  Figure  4. 
In  the  right  panel  we  show  the  corresponding  images  of  PSFs.  These  first  results  are  very  encouraging. 
The  true  and  estimated  MTFs  overlap  and  structure  of  the  PSF  is  well  represented  (note  the  AO  correction 
zone,  roughly  14  pixels  in  size,  in  panel  c).  The  results  are  more  than  satisfactory  given  that  the  fiber  PSF 
was  taken  around  two  years  after  the  first  observations. 


Fig.  4.  Illustration  of  the  method’s  performance.  Left:  Black  line:  ground  truth,  i.e.  the  MTF 
corresponding  to  the  observations.  Green  diamonds:  estimated  MTF  of  the  turbulence  AO-residuals 
only.  Blue  triangles:  estimated  AO  MTF  (green  diamonds)  multiplied  by  the  perfect,  diffraction -limited 
MTF.  Red  squares:  AO  MTF  multiplied  by  the  MTF  containing  the  telescope  and  static  error  terms 
(taken  from  a  fiber  PSF).  Right:  (a)  ground  truth  -  in  this  case  real  PSF,  (b)  fiber  PSF  (no  turbulence, 
only  static  errors),  (c)  result  of  PSF  reconstruction  after  inclusion  of  the  perfect,  diffraction -limited 
MTF,  (d)  result  of  PSF  reconstruction  after  inclusion  of  the  information  from  the  fiber  PSF. 


Given  these  very  encouraging  results,  we  proceeded  to  tackle  a  much  more  challenging  but  also  much 
more  realistic  scenario  of  an  extended  object  which  changes  throughout  the  observations.  Data  was  kindly 
provided  by  AFRL  staff  involved  in  this  project.  The  set  comprises  100  short  exposures  (texp  =  15  ms)  of  a 
defunct  US  remote-sensing  satellite  Seasat  taken  with  the  3.5-m  telescope  at  the  Starfire  Optical  Range, 
New  Mexico.  Observations  were  carried  out  in  the  laser-guide  star  AO  mode,  with  24x24  actuators  on  the 
deformable  mirror.  Wavelength  of  the  observations  was  0.8  pm.  It  is  estimated  that  the  Strehl  ratio  of  the 
observations  is  in  the  range  10-30  %.  Figure  5  shows  some  frames  from  the  sequence. 
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Fig.  5.  Example  frames  from  the  Seasat  sequence.  Note  the  appearance  and  disappearance  of  glints, 
change  in  aspect  and  local  illumination  of  the  satellite  and  the  rotation  with  respect  to  the  observer. 


There  are  two  challenging  aspects  associated  with  the  application  of  the  described  PSF  reconstruction 
method  on  this  data.  Firstly,  the  object  is  rapidly  changing:  the  satellite  is  rotating,  tumbling,  changing 
aspect,  glinting,  etc.  Secondly,  in  contrast  to  celestial  objects,  the  satellite’s  overhead  movement  is  very 
fast  and  the  level  of  AO  correction  changes  throughout  the  sequence.  In  order  to  obtain  meaningful  PSFs, 
the  sequence  was  therefore  divided  into  subsequences.  This  was  a  trade-off  study:  for  shorter 
subsequences  the  object  is  more  constant  but  PSF  estimation  is  less  accurate  because  of  small  samples. 
Longer  subsequences  have  better  signal-to-noise  ratio  on  the  estimators  but  the  variability  of  the  object 
leaks  through  Equations  (3)  and  (4).  After  a  trade-off  study  a  subsequence  length  of  25  frames  was  found 
optimal  for  these  observations.  Figure  6  shows  input  image  centered  in  the  middle  of  one  such 
subsequence,  the  reconstructed  PSF,  as  well  as  the  results  of  various  deconvolution  algorithms  which  were 
supplied  with  this  PSF.  The  algorithms  did  not  update  the  PSF;  they  simply  used  it  as  it  was  (non-blind 
mode).  The  improvement  in  image  quality  with  regard  to  the  original  data  shows  that  the  reconstructed 
PSF  is  useful.  Naturally,  a  more  advanced  deconvolution  algorithm  could  update  this  first-guess  PSF 
(“myopic”  mode)  and  offset  the  errors  which  arise  in  the  PSF  reconstruction  part  due  to  e.g.  changes  in  the 
object  or  estimation  noise. 


input  data 
sum  of  3  frames 

✓ 

deconvolved,  R-L 

* 

S' 

deconvolved,  AWMLE 

S' 

deconvolved,  TV 

* 

S' 

PSF  obtained  from  25  frames 

* 

* 

w 

* 

p 

Fig.  6.  Left:  input  data  whereby  three  short-exposure  frames  were  summed  to  increase  the  signal-to- 
noise  ratio.  Second  from  left:  reconstructed  PSF.  Center:  input  data  deconvolved  with  the  reconstructed 
PSF  and  the  Richardson-Lucy  algorithm  (R-L).  Second  from  right:  deconvolution  result  with  the 
AWMLE  algorithm  [11].  Far  right:  deconvolution  with  the  total  variation  regularization. 
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How  can  this  new  method  be  used  to  tackle  the  problem  of  anisoplanatic  deconvolution?  In  the  third  year 
of  the  project  it  will  be  combined  with  the  image  decomposition  approach  developed  in  the  first  year:  we 
studied  ways  of  optimally  describing  an  anisoplanatic  image  in  terms  of  the  principal  components  of  the 
spatially-varying  PSF.  Now  we  have  a  way  to  find  the  local  PSFs.  Rewriting  Equation  (1)  for  the 
anisoplanatic  case  we  have: 


i(x)  =  o(v)®h(v,3c)  +  n(3c)  (14) 

This  time  we  have  deliberately  included  noise  n(x)  in  the  equation.  The  symbols  x  and  v  are  2-D  vectors 
denoting  image-,  and  object-space  coordinates,  respectively.  Note  that  the  PSF  is  assumed  to  change  with 
v  and  x.  We  assume  that  the  set  of  reconstructed  PSFs  can  be  decomposed  into  J  principal  components 
h;  via  singular  value  decomposition: 


h(v,*)«Zay.(v)-h/(x) 


j= i 


(15) 


Then  Equation  (14)  becomes: 


i(V>«X[(«, -o)®h7.](x)  +  n(x)  (16) 

1=1 

The  deconvolution  process  would  then  involve  inverting  Equation  (16),  possibly  under  some 
regularization,  in  order  to  obtain  the  object  o(x) .  Because  in  practice  the  PSFs  are  most  probably  only 

going  to  be  recoverable  at  a  limited  number  of  locations,  the  coefficients  are  going  to  be  sparse.  Some 

interpolation  would  have  to  be  employed  in  order  to  obtain  a  map  of  covering  the  whole  object.  A 

regularized  inversion  of  Equation  (16)  would  then  allow  to  offset  the  three  sources  of  error:  PSF 
reconstruction  error,  coefficient  spatial  interpolation  error  and  the  influence  of  noise. 

4.2  Marginal  blind  deconvolution  applied  to  simulated  data 

For  testing  of  the  marginal  criterion  on  images  of  satellites,  a  set  of  6  AO  long-exposure  PSFs  were 
simulated  analytically  using  the  PAOLA  package  [12].  The  simulations  were  carried  out  for  the  3.5-m 
Starfire  Optical  Range  telescope  located  at  a  site  of  rather  bad  seeing  in  New  Mexico.  A  common  value  of 

1.2  "  for  seeing  was  adopted  (the  seeing  angle  is  defined  as  the  full-width-at-half-maximum  of  the  long- 
exposure  PSF  taken  through  turbulence  at  500  nm  wavelength  through  an  arbitrarily  large  telescope).  With 
Nyquist-sampled  PSFs  (two  pixels  per  A/D)  a  256  x  256  pixel  field  of  view  (FOV)  covers  7.5".  The 
magnitude  of  the  guide  star  (GS)  was  set  to  10  in  the  visible,  the  overall  transmission  from  GS  to  detector 
was  set  to  15%.  The  zenith  angle  was  set  to  30  °.  The  wavelength  of  the  observations  was  set  to  1  pm.  The 
AO  system  has  many  parameters  but  most  of  them  are  outside  the  scope  and  interest  of  this  report.  The 
most  important  parameter  is  the  density  of  actuators  on  the  deformable  mirror  and  simulations  followed 
the  current  design  of  the  SOR  AO  system  with  24  X  24  actuators.  The  other  parameters  were  either  taken 
from  SOR  AO  design  documents  [13]  or  optimized  to  obtain  the  highest  Strehl  ratio  on  axis.  Moreover, 
simulated  PSFs  were  obtained  with  two  real  turbulence  profiles  which  had  been  obtained  in  the  course  of 
site  testing  in  the  context  of  the  European  Extremely  Large  Telescope  project  [14].  One  profile  has  strong 
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ground-level  turbulence  and  produces  little  anisoplanatism  while  the  contrary  is  true  for  the  second 
profile. 


Fig.  7.  Simulated  set  of  PSFs  in  anisoplanatic  conditions  for  two  different  atmospheric  profiles  (rows) 
and  three  different  positions  (columns) 

In  the  simulations,  firstly  a  GS  image  was  generated  and  put  in  the  top  right  comer  of  the  FOV.  Then,  a 
loop  over  positions  generated  a  PSF  for  the  desired  locations  in  the  FOV  according  to  its  field  angle  and 
orientation  with  respect  to  the  GS.  These  PSFs  were  then  stored  for  subsequent  convolution/deconvolution 
and  are  shown  in  Figure  7.  As  test  object  we  used  the  schematic  image  of  the  Operationally  Responsive 
Space-1  satellite  (ORS)  depicted  in  Figure  9.  A  test  image  was  produced  by  convolving  the  test  object 
with  one  of  the  six  PSFs.  Subsequently,  shot  noise  and  read  out  noise  (le'  RMS)  were  added  to  produce 
the  final  images.  Two  different  levels  of  noise  were  considered  assuming  different  mean  fluxes  for  the 
object,  equal  to  103  and  104  photons  per  pixel. 

Finally,  the  marginal  criterion,  with  and  without  positivity  on  the  object,  was  tested  over  these  images  to 
obtain  the  coefficients  that  synthesized  the  PSF,  to  assess  how  well  we  can  recover  the  true  PSF  out  of  the 
set  of  six  simulated  ones.  The  final  goal  is  to  recover  a  good  PSF  that  can  be  used  with  whatever  non-blind 
deconvolution  algorithm  we  wish  to  use.  Here,  we  have  used  the  natural  output  of  each  approach,  i.e.,  for 
the  marginal  criterion  with  positivity  an  explicit  object  is  obtained  at  each  iteration,  hence,  in  Figure  9  we 
show  the  object  at  the  last  iteration.  For  the  marginal  criterion  with  no  positivity  the  object  is  computed 
implicitly  during  the  minimization  process  by  means  of  the  Wiener  solution,  i.e.,  there  is  no  object  at  the 
last  iteration.  Therefore,  in  order  to  produce  the  final  result  we  used  the  recovered  PSF  and  object  PSD 
together  with  a  Wiener  filter. 

Table  1  and  Figures  8  and  9  show  the  results  of  the  marginal  estimator  when  the  image  was  simulated  for 
the  PSF  number  3  (top  row,  third  column)  in  Figure  7.  This  is  a  representative  result  for  the  rest  of  the 
tested  PSFs.  Table  1  shows  that  the  marginal  criterion  is  a  powerful  tool  to  estimate  the  correct  PSF  for 
long  exposure  AO  observations,  from  a  single  image.  The  coefficient  corresponding  to  the  PSF  which 
created  the  image  was  correctly  identified  by  the  estimator,  with  more  robustness  when  the  positivity  on 
the  object  was  applied.  Figure  8  shows  the  recovered  PSFs  and  Figure  9  the  final  reconstructed  object  with 
the  corresponding  PSF,  either  with  or  without  positivity. 
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Table  1.  Results  of  the  marginal  estimator  applied  to  test  data  created  with  the  PSF  in  the  top  row  and  third  column  in  Figure  7. 
The  input  data  is  shown  in  Figure  9 
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Fig.  8.  Top  row:  PSF  estimation  with  the  marginal  criterion.  Bottom  row:  PSF  estimation  with  the 
marginal  criterion  and  positivity  on  the  object.  Leftmost  panels:  ground-truth  PSFs.  Middle  panels: 
estimated  PSFs.  Rightmost  panels:  difference  between  the  real  PSF  and  the  estimated  one  (log  scale). 
Object  mean  flux  equals  103  photons  per  pixel. 
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Fig.  9.  Top  row:  results  at  high  SNR,  object  mean  flux  equal  to  104  photons/pixel.  Bottom  row:  results  at 
low  SNR,  object  mean  flux  equal  to  103  photons/pixel.  First  column:  ground  truth  object.  Second 
column:  simulated  image  after  convolution  with  AO  long-exposure  PSF  in  Figure  7  (top  row,  third 
column).  Photon  noise  plus  read  out  noise  (le‘  RMS)  were  consecutively  added.  Third  column:  results 
from  the  marginal  criterion  with  positivity  on  the  object.  An  explicit  object  is  obtained  at  each  iteration; 
these  results  show  the  restored  object  at  the  last  iteration.  Fourth  column:  results  from  the  marginal 
criterion  with  no  positivity  on  the  object.  An  implicit  object  is  computed  using  the  Wiener  solution  for  a 
given  PSF,  hence,  the  final  object  is  obtained  using  the  final  PSF  and  the  Wiener  filter. 

5.  Conclusions 

In  the  second  year  of  the  project  we  have  introduced  a  novel  way  of  extracting  a  PSF  from  sequence  of 
images  of  arbitrary  objects.  The  method  is  applicable  to  astronomical  imaging  with  AO  and  to  imaging 
with  small  apertures  and/or  long  wavelengths.  Especially  interesting  will  be  the  application  of  the  new 
approach  to  estimation  of  spatially-varying  PSF  in  single-conjugate  AO  imaging.  Based  on  experience 
with  real  data  from  Starfire  Optical  Range  we  posit  that  the  method  will  work  with  less  than  100  images. 
There  is  also  nothing  preventing  the  use  of  the  method  with  long  exposures:  the  only  thing  that  will 
change  will  be  then  the  number  of  random  phasors  N  which  will  correspondingly  increase.  Practical  issue 
might  arise  in  that  both  contrast  and  skewness  diminish  towards  zero  when  A  increases  (see  Equation  (6)). 
This  might  lead  to  degeneracy  of  the  solution  and  theoretical  work  is  already  underway  to  find  more 
accurate  models  for  speckle  statistics  [15].  This  and  other  issues,  especially  pertaining  to  anisoplanatic 
deconvolution  based  on  singular  value  decomposition  of  the  reconstructed  local  PSFs,  will  be  explored  in 
Year  3  of  the  project. 

We  also  have  extended  the  marginal  blind  deconvolution  from  the  context  of  retinal  imaging  towards  the 
identification  of  satellites  observed  with  single-conjugate  AO  systems.  The  marginal  criterion  has  been 
modified  to  include  constraints  on  the  object,  such  as  positivity,  in  order  to  make  the  PSF  estimation  more 
robust  in  the  presence  of  noise  [16].  This  approach  can  be  extended  to  the  anisoplanatic  case  by  its 
application  to  different  positions  in  the  FOV  as  an  essential  building  block  of  a  global,  shift-variant  image 
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restoration  method.  Additionally,  we  plan  to  create  a  large  dictionary  of  PSFs,  describing  different 
atmospheric  conditions,  from  which  marginal  blind  deconvolution  will  extract  the  correct  linear 
combination  of  PSF  coefficients  at  each  position.  In  order  to  help  the  convergence  of  the  algorithm 
different  constraints  over  the  PSF  can  be  applied,  such  as  sparsity  on  the  coefficients  through  an  11 
regularization  in  order  to  force  the  real  PSF  to  a  be  a  linear  combination  of  only  a  few  of  the  PSFs  within 
the  dictionary,  and/or  constraining  the  shape  and  orientation  of  the  PSF  to  correspond  with  its  position  in 
the  FOV. 

Main  research  thrust  in  Year  3  will  be  devoted  to  development  of  a  novel  kind  of  simulation  tool  which 
will  allow  for  rapid  generation  of  anisoplanatic  image  sequences  for  three  scenarios:  uncompensated, 
horizontal-path  imaging  with  very  small  isoplanatic  angles,  single-conjugate  AO-compensated  vertical 
imaging,  and  air-to-ground,  “looking-down”,  slant-path  scenario.  Some  work  on  the  theory  behind  such  a 
simulation  has  been  carried  out  in  Year  2  of  the  project  [17]. 

The  methods  of  PSF  estimation  and  anisoplanatic  deconvolution  developed  in  the  course  of  the  project 
until  now  would  then  be  tested  on  these  simulated  sequences. 
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7.  List  of  Symbols,  Abbreviations,  and  Acronyms 

AO  -  adaptive  optics 

AMOS  -  Air  Force  Maui  Optical  and  Supercomputing 
AWMLE  -  Adaptive  Wavelet  Maximum  Likelihood  Estimator 
FOV  -  field  of  view 
GS  -  guide  star 

MTF  -  modulation  transfer  function 

ORS  -  Operationally  Responsive  Space- 1  satellite 

OTF  -  optical  transfer  function 

PSF  -  point-spread  function 

R-L  -  Richardson-Lucy  algorithm 

RMS  -  root  mean  square  (standard  deviation) 

SNR  -  signal-to-noise  ratio 
TV  -  total  variation  constraint 
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